library(foreign)
library(ggeffects)



setwd("/Users/joelsievert/Dropbox/Research/Legislative Responsiveness/National Road/data/PRQ Replication")


####################################
##votes to transfer road to states##
###################################
h21 <- read.csv("h21_extend.csv")
h23 <- read.csv("h23_extend.csv")

#national road

m1 <- glm(rc_natl ~ jackson  + dist_road + dist_natl + slave_p + distrib_cmte, family = binomial, data = h21)

#maysville
m3 <- glm(rc_mays ~ jackson + dist_road + dist_mays + slave_p + distrib_cmte, family = binomial, data = h21)


##23rd
m5 <- glm(rc ~ jackson   + dist_road + dist_natl + slave_p + distrib_cmte , family = binomial, data = h23)


p1 <- ggpredict(m1, terms = "slave_p[0, 0.04, 0.08, 0.12, 0.16, 0.2, 0.24, 0.28]", condition = c(jackson = 1, dist_road = 2.1, dist_natl = 3.3, distrib_cmte = 0), ci_level = 0.9) #, vcov_type = sandwich:vcovHC)
p2 <- ggpredict(m3, terms = "slave_p[0, 0.04, 0.08, 0.12, 0.16, 0.2, 0.24, 0.28]", condition = c(jackson = 1, dist_road = 2.1, dist_mays = 5.3, distrib_cmte = 0), ci_level = 0.9) #, vcov.type = sandwich:vcovHC)
p3 <- ggpredict(m5, terms = "slave_p[0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.15, 0.175, 0.20, 0.225, 0.25]", condition = c(jackson = 1, dist_road = 2.38, dist_natl = 9.6, distrib_cmte = 0), ci_level = 0.9) #, vcov.type = sandwich:vcovHC)

#p4 <- ggpredict(m3, terms = "slave_p[0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.15, 0.175, 0.20, 0.225, 0.25]", condition = c(jackson = 0, dist_road = 2.38, dist_natl = 9.6, distrib_cmte = 0), ci.lvl = 0.9, vcov.type = sandwich:vcovHC)


par(mfrow = c(1,3))

#Washington Road
par(mar = c(5, 4, 3, 2), bg = "white")
plot(p1$x, p1$predicted, ylim = c(0, 1), xlim = c(0, 0.28), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(0, 0.28, .01), col = 'grey85', lwd = 80)
abline(v = seq(0, 0.28, 0.04), h = seq(0, 1, 0.1), col = 'white')
lines(p1$x, p1$predicted, lwd = 3)
lines(p1$x, p1$conf.low, lwd = 3, lty = 2)
lines(p1$x, p1$conf.high, lwd = 3, lty = 2)

#code puts axis text and an angle
axis(side = 2, at = seq(0, 1, 0.1), labels = seq(0, 1, 0.1), las = 2, lwd = 0 ) 
axis(side = 1, at = seq(0, 0.28, 0.04), lwd = 0)

mtext("Slave Population (%)", side = 1, line = 2.75)
mtext("Predicted Probability of Support", side = 2, line = 2.5)

mtext("Washington Road", side = 3, line = 0.5, cex = 1, font = 2)


#Maysville Road
par(mar = c(5, 4, 3, 2), bg = "white")
plot(p2$x, p2$predicted, ylim = c(0, 1), xlim = c(0, 0.28), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(0, 0.28, .01), col = 'grey85', lwd = 80)
abline(v = seq(0, 0.28, 0.04), h = seq(0, 1, 0.1), col = 'white')
lines(p2$x, p2$predicted, lwd = 3)
lines(p2$x, p2$conf.low, lwd = 3, lty = 2)
lines(p2$x, p2$conf.high, lwd = 3, lty = 2)

#code puts axis text and an angle
axis(side = 2, at = seq(0, 1, 0.1), labels = seq(0, 1, 0.1), las = 2, lwd = 0 ) 
axis(side = 1, at = seq(0, 0.28, 0.04), lwd = 0)

mtext("Slave Population (%)", side = 1, line = 2.75)
#mtext("Predicted Probability", side = 2, line = 2.75)

mtext("Maysville Road", side = 3, line = 0.5, cex = 1, font = 2)


#New Orleans Route
par(mar = c(5, 4, 3, 2), bg = "white")
plot(p3$x, p3$predicted, ylim = c(0, 1), xlim = c(0, 0.25), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(0, 0.25, .01), col = 'grey85', lwd = 80)
abline(v = seq(0,0.25, .05), h = seq(0, 1, 0.1), col = 'white')
lines(p3$x, p3$predicted, lwd = 3)
lines(p3$x, p3$conf.low, lwd = 3, lty = 2)
lines(p3$x, p3$conf.high, lwd = 3, lty = 2)

#code puts axis text and an angle
axis(side = 2, at = seq(0, 1, 0.1), labels = seq(0, 1, 0.1), las = 2, lwd = 0 ) 
axis(side = 1, at = seq(0, 0.25, 0.05), lwd = 0)

mtext("Slave Population (%)", side = 1, line = 2.75)
#mtext("Predicted Probability", side = 2, line = 2.75)

mtext("New Orleans Route", side = 3, line = 0.5, cex = 1, font = 2)



